TP Simulations de Variables Aléatoires

import numpy as np
from matplotlib import pyplot as plt
import numpy.random as nprd
from scipy.stats import norm,cauchy
from IPython.display import clear_output
normal_d,cauchy_d = norm(),cauchy()
from time import time

Partie 1

A partir de l'inverse généralisée de la fonction de répartition

L'inverse généralisé de F est donné par : $F^-1 = -\dfrac{1}{\lambda} ln(1-u)$

def simulation_echantillons(n=1000,Lambda=2,norm=1):
    tab_echantillons=[]
    for i in range(n):
        u=np.random.uniform()
        tab_echantillons.append(((-1/Lambda)*np.log(1-u))/norm)
    return tab_echantillons

plt.style.use('ggplot')
particules=simulation_echantillons(1000,2)
sum=np.sum(particules)
particules=[p/sum for p in particules]
plt.hist(particules,bins=100)
plt.title("Echantillons d'une variable aléatoire exponetielle")      

Text(0.5, 1.0, "Echantillons d'une variable aléatoire exponetielle")

norm = lambda a : a/1000
x=np.linspace(0,4,1000)
f=2*np.exp(-2*x)
plt.figure()
plt.plot(x,f,label='densité de probabilité exponentielle pour lambda = 2')
particules=simulation_echantillons(1000,2)
plt.hist(particules,bins=100,density=True)
plt.ylim(0,2)
plt.legend()

<matplotlib.legend.Legend at 0x2d571334d88>

On constate que les échantillons normalisés approchent l'aire sous la courbe de densité. Pour N tend vers l'infini on montre que la somme des aires de l'histogramme tend vers l'aire sous la courbe.
Le problème de cette méthode est qu'elle demande une grande quantité d'achantillons pour estimer l'aire sous la courbe.

Partie 2

M=5
N_sample = 10000
def Accept_Reject(M=5,N_sample = 10000):
    X_gen=[]
    for j in range(N_sample):
        u = nprd.uniform()
        x = nprd.standard_cauchy()
        while 1/M*normal_d.pdf(x)/cauchy_d.pdf(x) < u:
            u=nprd.uniform()
            x=nprd.standard_cauchy()
        X_gen.append(x)
    return(X_gen)



    

X_gen = Accept_Reject()
plt.style.use("ggplot")
plt.hist(X_gen,bins= N_sample//200)
plt.title(" Méthode Accept - Reject")
plt.show()

t=np.linspace(normal_d.ppf(0.0001),normal_d.ppf(0.9999),10000)
X=normal_d.pdf(t)
plt.plot(t,X)
plt.title(" Loi normale centrée réduite")

Text(0.5, 1.0, ' Loi normale centrée réduite')

plt.figure(figsize=[2*5.4, 1*2.8])
plt.hist(X_gen,bins= N_sample//200,color="grey",label="Sample Data",density=True)
plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy ")
plt.legend()
plt.show()

On constate qu'on approche bien une loi Normale centrée réduite. Ainsi, la méthode $\textit{Accept-reject}$ est une bonne méthode pour générer de nombres aléatoires d'une densité connue à l'aide d'un autre générateur de nombres aléatoires, aussi de densité connue.
On pourra discuter de l'influence de $M$ qui correspond à l'inverse de la probabilité que x soit accepté est égale.

plt.figure()
X_norm=normal_d.rvs(size=N_sample)
plt.hist(X_gen,bins= N_sample//100,color="blue",label="Normal With Cauchy",alpha = 0.4,density=True)
plt.hist(X_norm,bins= N_sample//100,color="red",label="Normal",alpha = 0.3,density=True)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy et \n génération d'une loi normale centrée réduite")
plt.legend()
plt.show()

Qualitativement, on remarque bien que les deux histogrammes se supperposent et on en déduit qu'on a bien un générateur d'une loi normale à partir d'un générateur de nombres aléatoires suivant une loi de Cauchy

Discussion selon le paramètre $M$

m = [1,2,4,8,16]
plt.figure()
T_exec=[]
for _,M in enumerate(m):
    plt.figure(figsize=[2*5.4, 1.5*2.8])
    t_deb=time()
    X_gen = Accept_Reject(M)
    t_exec = time()-t_deb
    plt.hist(X_gen,bins= N_sample//20,color="grey",label="Sample Data",density=True)
    
    T_exec.append(t_exec)
    print(f"Temps d'exécution pour M = {M} vaut {t_exec} secondes" )
    plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
    plt.title(f" Méthode Accept - Reject avec $M={M}$")
    plt.show()
plt.figure()
plt.plot(m,T_exec,color = "black")
plt.xlabel("$M$");plt.ylabel("Temps d'éxécution");plt.title("Temps d'éxécution en fonction de $M$")
plt.semilogx()
plt.show()

Temps d'exécution pour M = 1 vaut 4.822746276855469 secondes

<Figure size 432x288 with 0 Axes>

Temps d'exécution pour M = 2 vaut 7.48268723487854 secondes

Temps d'exécution pour M = 4 vaut 14.612187623977661 secondes

Temps d'exécution pour M = 8 vaut 30.970727682113647 secondes

Temps d'exécution pour M = 16 vaut 62.628419160842896 secondes

Plus M est grand, mieux on approxime la loi mais plus le temps d'exécution augmente

3 Algorithme de Box-Muller

Dans l'exercice, au lieu de regarder les résulats à 1 dimension, on utilisera plotly pour regarder les résulats en 2 dimensions.

N_sample = 10000
def Box_Muller(N_sample = 10000):
    X_gen,Y_gen=[],[]
    for j in range(N_sample):
        u1 = nprd.uniform()
        u2=nprd.uniform()
        R= -2*np.log(u1)
        V=2*np.pi*u2
        X_gen.append(np.sqrt(R)*np.cos(V))
        Y_gen.append(np.sqrt(R)*np.sin(V))
    return(X_gen,Y_gen)

import seaborn as sns 
X,Y = Box_Muller(10000)

plt.figure()
plt.scatter(X,Y)
plt.show()

plt.figure()
sns.kdeplot(x=X, y=Y,cmap="Blues",fill=True, thresh=0 )
plt.show()

Pour utiliser plotly, installez le à l'aide de $\textit{pip install plotly}$.
Si vous avez besoin de sauvegarder une image, utilisez fig.write_image("Image.png"), mais installez avant kaleido avec $\textit{pip install -U kaleido}$.

import plotly.graph_objects as go
h,x1,y1,image = plt.hist2d(X,Y,bins=20,density=True)
fig = go.Figure(data=[go.Surface(z=h, x=x1, y=y1)])
# fig.write_image("Image.png")
fig.show()

plt.imshow(plt.imread("Image.png"))

<matplotlib.image.AxesImage at 0x2d516f8dd88>

On remarque qu'on génère bien une gaussienne en 2D à l'aide de la méthode de Box-Muller. Néanmoins, c'est une méthode plus couteuse en temps.

from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([0,0],np.identity(2))
x,y=np.mgrid[-3:3:.1, -3:3:.1]
z= np.dstack((x, y))

fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),go.Surface(z=h, x=x1, y=y1, opacity=0.4,colorscale='Reds',showscale=False)])
fig.write_image("Image3.png")
fig.show()

plt.imshow(plt.imread("Image2.png"))

<matplotlib.image.AxesImage at 0x2d516cf5f08>

plt.imshow(plt.imread("Image3.png"))

<matplotlib.image.AxesImage at 0x2d516d0c8c8>

La gaussienne 2D se confond avec l'histogramme 2D généré avec la méthode de Box-Muller. En relançant plusieurs fois le programme, on remarque qu'on conserve bien la coincidence entre l'histogramme et la surface de densité de la gaussienne bivariée.

Désormais, générons une loi normal de moyenne -3 et d'écart type 3.

def Box_Muller(N_sample = 10000,mu=-3,sigma=3):
    X_gen,Y_gen=[],[]
    for j in range(N_sample):
        u1 = nprd.uniform()
        u2=nprd.uniform()
        R= -2*np.log(u1)
        V=2*np.pi*u2
        X_gen.append(sigma*np.sqrt(R)*np.cos(V)+mu)
        Y_gen.append(sigma*np.sqrt(R)*np.sin(V)+mu)
    return(X_gen,Y_gen)

X,Y =  Box_Muller(N_sample = 5000,mu=-3,sigma=3)

from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([-3,-3],3**2*np.identity(2))
x,y=np.mgrid[-10:4:.1, -10:4:.1]
z= np.dstack((x, y))
h,x2,y2,image = plt.hist2d(X,Y,bins=30,density=True)

fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),
                      go.Surface(z=h, x=x2, y=y2, opacity=0.4,colorscale='Reds',showscale=False)])
fig.show()

On obtient ce résulat en remplaçant $\sqrt{R}\sin{V}$ par $\sigma\sqrt{R}\sin{V}+\mu$

4 Algorithme de Generation d'un echantillon normal multivarie de dimension n

from scipy.stats import multivariate_normal

def gen_N_mu_sigma(mu,sigma,d=1000):

    n,X_gen=len(mu),[]

    norm_nd = multivariate_normal(np.zeros(n),np.identity(n))
    A = np.linalg.cholesky(sigma)
    for j in range(d):
        Z = norm_nd.rvs()
        X_gen.append(Z@A+mu)
    return(np.array(X_gen))

import seaborn as sns
X= gen_N_mu_sigma([10,1],[[0.4,0.2],[0.2,0.4]])
plt.figure(figsize=(10,10))
sns.kdeplot(x=X[:,0],y=X[:,1])
plt.show()

mu = np.array([0 ,50 ,100 , 50 ,100 ,200])
sigma = np.array([[11, 10, 5, 9, 4, 2]
        ,[10 ,13, 9, 15, 5, 3]
        ,[5, 9 ,15, 11, 3, 1 ]
        ,[9, 15 ,11 ,21, 6, 4]
        ,[4 ,5 ,3 ,6, 5, 1]
        ,[2 ,3 ,1 ,4, 1, 1]])


X= gen_N_mu_sigma(mu,sigma, 10000)

plt.figure(figsize=(10,5))
sns.kdeplot(X)
plt.show()

for i in range(6):
    print(f"La variance de Z{i} est de {np.diag(sigma)[i]}")
    

La variance de Z0 est de 11
La variance de Z1 est de 13
La variance de Z2 est de 15
La variance de Z3 est de 21
La variance de Z4 est de 5
La variance de Z5 est de 1

Plus la variance est grande, plus la courbe de densité est étalée. C'est ce qu'on remarque sur le graphique ci dessus.

L'amplitude max est donnée par $$\dfrac{1}{\sqrt{2\pi}\sigma}$$
Ainsi, plus sigma est petit, plus l'amplitude max est grande.

import pandas as pd


sns.pairplot(pd.DataFrame(X,columns=[f"Z{i}" for i in range(6)]))
plt.show()

Indice, Corr = [],[]
for i in range(1,len(sigma)):
    for j in range(i):
        Corr.append(sigma[i,j]/(np.sqrt(sigma[i,i])*np.sqrt(sigma[j,j]))) 
        Indice.append(f"{i},{j}")
pd.DataFrame(np.array([np.choose(np.flip(np.argsort(Corr)), np.array(Indice)),np.flip(np.sort(Corr))]).T,columns=["Indices","Corréalation"]).head(5)

	Indices	Corréalation
0	3,1	0.9078412990032038
1	5,3	0.8728715609439696
2	1,0	0.8362420100070909
3	5,1	0.8320502943378437
4	2,1	0.6445033866354896

On peut voir les différentes dispersions selon les diverses sous espaces.
Les variable $Z_3$ et $Z_1$ sont fortement corrélée et on peut le remarquer sur le graphique pairplot ci dessus. On a pu calculer les coefficients de corrélation : $$\rho_{i,j} = \dfrac{\Sigma_{i,j}}{\sigma{i}\sigma{j}}$$

5 Echantillonner suivant une loi de Bernouilli

def Bern(p,d=1000):
    B_gen=[]
    for j in range(d):
        u = nprd.uniform()
        B_gen.append(0 if u <p else 1)
    return(np.array(B_gen))

B = Bern(0.7)
plt.hist(B,density=True)
plt.show()

nb_zeros = (len(B)-np.sum(B))/len(B)
print(f"La fréquence du nombre de zéros est de {nb_zeros}")

La fréquence du nombre de zéros est de 0.696